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Abstract 


Holographic Particle Image Velocimetry (HPIV) is a novel technique for measuring 
the complete fluid flow around a body. Advances in computing power make this 
technique practicable for the first time. 

Currently popular techniques for experinentally determining fluid flow around 
a test. body rely on measuring the flow at a single point and moving the sample 
point during the experiment. This implicitly time averages the data. HPIV can 
measure the flow field in a volume nearly instantaneously by sequentially “holo- 
graphing” the seeded flow field. 

Until recently the only practical method of processing the HPIV data was to 
physically reconstruct the optical image from the hologram and scan through the 
resulting volume using microscope optics. This is a very slow process. 

Recent. advances in computing power make it possible to directly analyze the 
digitized holograph (digitized using a high resolution scanner) This reduces the 
post-processing time from several hours to less than 1-hour with much less expen- 
sive equipment. 

The primary contribution of this thesis is an advanced multi-platforin graphical 
user interface hologram synthesizer for use in calibrating digital HPIV reconstruc- 
tion algorithms. 
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Chapter 1 


Introduction 


1.1 Overview 


This thesis presents an advanced multiplatform graphical user interface hologram 
synthesizer to produce calibration holograms of seeded flow fields. The synthe- 
sizer, HoloSynth, runs under X11R6 with OpenGL, and Microsoft Windows NT. 
Complete source code is contained in Appendix A. A brief introduction to parti- 
cle holography and holographic particle velocimetry is presented as well as results 


from computational reconstruction of the hologram generated by the synthesizer. 


1.2 Motivation 


Until very recently, particle image velocimetry (PIV) was limited to measuring 
fluid flow velocities only in a single plane. A sheet of laser light was passed 
through the flow which contained neutrally buoyant particles. A double exposure 
of the particles onto either film or a Charge Coupled Device (CCD) camera allows 
the two velocity components in the illuminated plane to be determined. 

The PIV method has recently been extended to three dimensions by Katz et. al. 
1] by making a hologram of the flow field. Three dimensional Holographic Particle 
Image Velocimetry (HPIV) uses either in-line or off-axis holographic techniques 
(see Section 2.2). Once the hologram has been recorded, the image is reconstructed 


by reversing the recording process and using a three dimensional optical scanner to 
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digitize the image for analysis by computer. This process can take several hours. 

An alternative to scanning the reconstructed image is to digitize the hologram 
and allow a computer to reconstruct the image m software. The recent advances i 
Computmg power should make this process far faster than the mechanical /optical 
scanning process now used. 

Digitizing the hologram can be done in two ways. ‘The mitial method avoids 
the three dimensional scanning of the reconstructed image by using a high perfor- 
urance flat bed scanner to digitize au enlargement of the transparency produced 
in the normal holographic process. Alternately, a high resolution scanner can be 
used directly on the hologram transparency. The flat bed scanning process is well 
understood, since the transparency can be enlarged to overcome resolution prob- 
lems witl: the scanner. However, it still requires recording the hologram on filin 
with the attendant film processing. 

A more advanced technique would be to capture the hologram using a Charge 
Coupled Device (CCD) Camera. The hologram would be immediately available 
ui digital form, and no film or processing would be mvolved. 

Regardless of the hologram digitization technique the computational machinery 
for reconstructing the image must be calibrated. It far easier to produce a synthetic 
hologram for calibration than to set up a tightly controlled experiment and actually 


record a hologram. 


1.3. Requirements 


The initial goals of computationally reconstructed HPIV are to be able to accu- 
rately determine the flow field in a cubic volume a few inches on a side. The fluid 
speeds will be typical of those used in hydrodynamics research. The quantifiable 


requirements are: 


e Must be able to resolve a 10yni particle. This corresponds to the niost 
conveniently sized neutrally buoyant seed particles used for Laser Doppler 


Anemoinetry. 
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e Accurate coverage within a 125cm* (5cim cube) volume. This will rise as the 


technique is perfected, and computing power increases. 


e Unambiguously detect fluid speeds < 10m/s. This accounts for free stream 


plus allowance for turbulence and model scale factors. 


The first and second requirements dictate the size and resolution of the recording 
device. 

Since holography is essentially a linear process, multiple exposures can be used 
instead of a rapid succession of exposures. So rather than design a system to 
rapidly capture two separate holograims. the recording medium is exposed two 
or more times without being altered between exposures. This results in multiple 
real images for each particle in the volume, the analysis process must identify the 
multiple unages with a particular particle. ‘This alleviates the need for a high 
speed film traversal mechanism and removes CCD output data rate as a concern. 

The third requirement sets the pulse lengths and pulse intervals of the laser 
used to generate the hologram. If the pulses are too long the diffraction patterns 
will be smeared from the particle motion. The greater the interval between pulses, 
the harder it will be to uniquely identify particle images with a single particle. The 
pulse length and interval requirements determine the recording device sensitivity 


and laser output power. 





12 


Chapter 2 


Review of Holographic Technique 


2.1 Overview 


A holograi is a two dimensional recording that can produce a three dimensional 
image. Typically electromagnetic (light) waves are recorded. but recording of 
acoustic waves and matter waves have also been made[2]. Electromagnetic holog- 
raphy differs from normal photography in that it encodes phase information as 
well as intensity information on the recording mediun. 

Recording media respond only to the intensity of the incoming wavefront, not 
its phase. To encode the phase using normal media, holography depends on using a 
coherent reference beam to interfere with a scattered or diffracted wave, or “object 
wave’, from the objects under observation. As a result of this imterference, the 
intensity at the recording plane depends on both the phase and amphtude of the 
object wave. 

The original image is reproduced sunply by illuminating the hologram with 
coherent hght using geometry similar to that used to record the hologram. If the 
wavelength used for reconstruction differs from the wavelength used for recording. 
the image will be enlarged or reduced. 

The reconstruction process, while simple, is not amenable to rapid quantitative 
data collection. Quantitative analysis generally involves using three dimensional 


optical scanning of the reconstructed nage. The ability to computationally re- 





2.2. IN-LINE (GABOR) HOLOGRAPHY 13 


T 
Scattered Wave 4 






Monochromatic 
Point Source 
Transmitted Wave Recording Plane 


Figure 2-1: Basic in-line hologram geometry 


construct the image would drainatically increase the utility of holographic mea- 


surement techmiques. 


2.2 In-line (Gabor) Holography 


In-line holography, or Gabor Holography after its inventor, describes a holographic 
apparatus in which the reference beam is usually oriented perpendicular to the 
recording media and illuminates the objects from behind. It 1s not a useful tech- 
nique for objects that are large with respect to the recording medium, but if one 
is interested primarily in the three dimensional position of small objects and not 
their three dimensional shape, in-line holography can be very effective. It is com- 
monly used to size particulates and aerosols suspended in fluids, and it has been 
used for traditionally reconstructed HPIV{3}. 

Figure 2-1 shows the basic geometry and coordinate system for recording a 
plane in-lme hologram. The object is illuminated by the reference source. The 
diffracted wave, or object wave, interferes with the wave from the reference source. 
or reference wave. This interference pattern is recorded by the medium in the 
recording plane and constitutes the hologram itself. 

If the reference and object complex wave amplitudes measured at the recording 


plane are given by R (zx, y) and O (az, y) respectively, the intensity at the recording 
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Figure 2-2: In-line hologram reconstruction geometry 


plane is given by, 
a = = 2 
AE ye — |R (a, y) +O (x, y)| pz aly) 


Since the recording plane is perpendicular to the plane reference wave, the 


complex amplitude of the reference wave is a constant, r. Expanding (2.1) yields 


1(z,9) = |r+O(2,y)[° 


r+ |O(a,y)| +rO (2, y) + rO* (2,9) (2.2) 


The intensity is recorded in any acceptable manner, usually as a positive trans- 
parency on film. Most reasonable recording media respond nearly linearly to in- 
cident wave intensity. The transmittance, t (2, y), of a positive transparency can 


be adequately described by 


ile) 
ty — BT |r? + |O(a,y)/? + rO(e,y)+rO"(F,9)] (23) 


t (x,y) 


where 7 1s a constant background transmittance, T is the exposure time, and (3 
is a constant of the recording medium. 
Figure 2-2 shows an in-line hologram being physically reconstructed. The 


wavetront transmitted by the hologram is siniply the incident wave multiplied by 
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the transmittance of the hologram and can be written, 


v(ey) = rt(ay) 
r(m% + 8Tr*) + BTr |O io y)| 
+8Tr°O (2, y) + BT r°O* (2, y) (2.4) 


=) 


I 


The first term in (2.4) represents the uniformly attenuated reference wave. The 
second term is normally very small, since typically O(z,y) <r. It is called the 
“halo” and under proper conditions can be neglected. ‘The third term is identical 
to, within a multiplicative constant, with the original scattered wave. It forms an 
luinage 2, from the hologram on the side opposite the observer. This is the virtual 
image. which is what the observer perceives. The fourth term represents a wave- 
front similar to the originally scattered wavefront, but with opposite curvature. It 
converges to form a real image z, from the hologram, between the hologram and 
the observer.[4] 

The observer then sees the virtual image, with a superimposed, out-of-focus 
real image along with the diminished reference beam. The two images instead of 


one and are the biggest drawbacks to in-line holography. 


2.3 Off-axis (Leith-Upatneiks) Holography 


A simple way to separate the out-of-focus real image from the virtual image is to 
provide a separate path for the reference beam and the object wave. Figure 2-3 
shows the basic geometry for recording an off-axis hologram. 

Because the reference beam is no longer perpendicular to the recording plane, 


it cannot be simplified to a constant, rather, it can be expressed as, 


R(x, y) = rexp(j2n€,2) 


aS 
ho 
Cy 

——” 


where €, = (sin @)/A, since only the phase of the reference beam varies across the 


recording plane. 
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Figure 2-3: Off-axis hologram geometry 


In a fashion similar to the derivation of (2.2) the intensity at the recording 


plane is found to be, 


I(z,y)=17°|O (x,y)| + 2r|O (x, y)| cos [27£,2 + (2, y)| (Zao) 





Figure 2-4 shows the geometry for an off-axis hologram reconstruction. As- 
sunung the recording medium is hnear, the transmitted wave complex amplhtude 


is found in a similar manner, 


Tyr exp(j27E,£) 
+BTr |O (x,y) ° exp(j27€,2) 
+(Tr’O (x,y) 


w(x, y) 


+STr°O (x,y) exp(j47€,2) (2.7) 


Once again, the first term in (2.7) is simply the attenuated reference beam. 
The second term is a “halo” term, similar to the neglected term in (2.7). The 
third term is identical to the original scattered wave and forms a virtual iniage 
in the same position as the original object. The fourth term is, again, the real 
nage. this time shifted from the optical axis by approximately twice the angle 
the reference beam makes with the recording plane. Hence the out of focus real 


lage does not interfere with the vitual image. 
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Figure 2-4: Off-axis hologram reconstruction geometry 


At first blush it would seem off-axis holography is the answer to the virtual 
lunage problem for HPIV. If purely physical reconstruction is to be used, then off- 
axis holography is eminently practicable. The fatal flaw comes when considering 
digitizing the hologram for use in computational reconstruction. The reference 
bean now introduces a rapidly fluctuating intensity term that quadruples the 
resolution requirement of the recording medium and digitization system. For this 


reason off-axis holography will not be discussed further in this thesis. 


2.4 Diffraction 


Holography is fundainentally about diffraction. Goodman|5] provides a detailed 


derivation of the Huygens-Fresnel principle for a wave from constant z 


wenn) = 5 [f ten aot , 


wo 
OO 
—— 


where, 79) = 27 + (rn — 2)* + (yn — y)*, and x, and y, are coordinates in the 
recording plane. 

This expression, while physically meaningful, is not easily evaluated. Two 
fundamental simplifications are common depending on the geometry of the optical 
problem being analyzed. The approxiniations are the Fresnel and the Fraunhofer 


approxunations, valid in the near and far fields, respectively. 
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2.4.1 The Fresnel Approximation 


A simpler and more useful expression of the Huygens-Fresnel integral can be found 
: ‘ ; : ‘ . , ) 9 
by simplifying ro; using the paraxial approximation. 27 >> 2° 4+ y’. 


If b < 1, then, 
1 In 


Using (2.9) with b = (gu=2)" 7 (m=v)*, 


Tro. = ze + (ay = ee 5 (Yn rag y)? 








man) (2.10) 


Substitution of (2.10) into (2.8) yields little simplification unless carefully cho- 


( 
nN 
— 

NO | 
wae 
Ny | 
So ee 

bt 
ae 
NO | 
La ee 
N 


sen terms are dropped. For the rg, term in the denominator of (2.8) it is considered 
acceptable to drop all terms but z. However, small changes in the ro, occurring 
in the numerator can cause significant changes in the value of the exponential 
and cannot be omitted. These large changes are generated in part by the mul- 
tiplication of ro; by &, usually a very large quantity, and the sensitivity of the 
exponential to phase changes in its argument. The resulting simplification yields 


the Fresnel-Kirchoff Diffraction Integral: 





| eine T ake 
Uae) Sa a // U(x, y) exp ‘5. (xp — x) + (yp, — y) dx dy2.11) 
Cc 2m 





gkz _k 9 aS 2 
oi elas (th+¥R) / U(x, yet2= (2 +9") e-I5E @nz+9ny) dp ay{ 2.12) 
JAK 

= OO 
Where (2.12) differs from (2.11) only through factormg the term exp E (re + yp) 
outside the integral. This approximation holds when the observer is in the near 


field, or Fresnel zone, of the diffracting aperture. Typical HPIV geometry will 
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place the recording medium in the far field of the individual particles but in the 


near field of the ensemble. 


2.4.2 The Fraunhofer Approximation 


Additional siniplifications to the Fresnel Diffraction Integral can be made if the 


systen: meets the Fraunhofer approximation. 


Ge Aye) eos 





: 9 
z > 5 Clay) 
Sunplifving (2.12) yields, 
ae eee vi oe 
W(Tn, Yn) = DN ertlient) ff U(x, ye 9 x27 tua do dy (2.14) 
jrz 
— Oc 


which is valid in the far field. The integral is recognizable as a two dimensional 


Fourier transform. 


2.4.3. Convolution interpretation of Fresnel-Kirchoff Inte- 


gral 


The Fresnel-Kirchoff Integral, 


can easily be viewed as a convolution with kernel, 


esl kz 


Ne 





ho 
N | >> 


h.(r,y) = exp E - (x + )| (2.16) 


yielding 


Oc 


Y2(Thy Ya) = // U(x, y)he(Lh — ©, yn — y)dx dy = U(tp, yn) @hz(z,y) (2.17) 
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This expression determines the field diffracted from a planar object, or an enseible 
of co-planar planar objects. If a large number of very small non-coplanar objects 


are in the field of view, (2.17) can be extended: 


Wie Yn) =a SS CAG. Yn) & pe Con Yn) (2. 18) 


Th 


where the sum is taken over nm planes within the volume. Each object plane. 
U,,(xz, y). contains at least one diffracting object and is located z, from the record- 


ing plane. The planes need not be evenly spaced. 


2.5 Laser Characteristics 


2.5.1 Fundamentals 


At. the atomic level, light is always generated by an atom transitioning from one 
energy state, py, to a lower energy state, &;. When this transition occurs the atom 
emits a photon with frequency v = (Eo — E,)/h, where fh is Planck’s constant. 
This eiission is a probabilistic occurrence and may occur spontaneously provided 
the atom has sufficient internal energy. ‘This 1s spontameous emission. 

A second form of eniission, first suggested by Einstein in 1917, is called stim- 
ulated emission. In the case of stimulated emission, a photon collides with an 
excited atom. The mecident photon must have exactly the same energy as the 
transition the atom will undergo. When this condition is met, the atom will make 
the transition, emitting a second photon. The emitted photon is m phase with 
the incident photon and travels in the same direction. 

Stimulated emission is fundamental to laser operation. A material is induced to 
have a majority of its atoms in an excited state (called a “population inversion” ). 
Since the majority of the atoms in the material are excited to the same state. 
they have identical transition energies. Initial spontaneous emission by a small 
number of the atoms will stimulate emission in others. The stimulated photons 
will in turn stimulate further emission in a run-away chain reaction. This in itself 


will not cause light amplification, but if the material is mounted within an optical 
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cavity (two parallel reflectors), the stimulated emission will preferentially orient 
itself along tlie axis of the cavity. Radiation not parallel to the axis will be lost 
to the surroundings and not contribute to the electromagnetic wave building up 
in the cavity. Since the stimulated emission is always m phase, the radiation 
traveling along the cavity axis constructively reinforces and the amplitude rapidly 
increases forming a longitudinal vibration mode within the cavity. If an end of 
the cavity is only partially reflective, then some of the coherent amplified signal 


will be released, this is the usable laser beam. 


2.5.2 Establishing the Population Inversion 


The process of establishing the population inversion ts called ‘pumping’. Pumping 
can be accomplished using flash lamps, usual in the case of ruby lasers, or electrical 
discharge, as in the case of most gas lasers. Photochemical reactions can also 
stimulate the population inversion but that technique is very rare.[6] 

The cavity material characteristics fundamentally determine what type of pump- 
ing will be appropriate. Continuous beam lasers must obviously use a sustainable 
pumping technique, such as tungsten filament incandescent lamps, or the afore- 
mentioned electrical discharge. 

The length of the pulse m non-switched pulse lasers is directly related to the 
length of the pulse length of the pump. Flash lamp pumped ruby lasers typically 
provide pulses on the order of several milliseconds long. The length of the pulse 
also varies from pulse to pulse as the condition of the flash lamp changes. 

The minimum time between lasers pulses in non-switched lasers is directly 
related to the requirements of the pump. Flash lamps typically must cool imme- 
diately after use. The amount of cooling is directly related to the power output 
of the lamp, and can be on the order of minutes. This limits the pulse repetition 
frequency of the unswitched laser. 

Switching can reduce pulse lengths to femtoseconds (1x107!° seconds) and 


pulse intervals similarly. Switching is discussed in some detail in Section 2.5.5. 
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2.5.3 Energy Levels 


When common laser materials are pumped, the resulting atomic population usu- 
ally has many more than one possible transition state available. The lkelihood 
that an atom undergoes a particular transition depends on the quantum mechan- 
ics of the particular material. In any event, when the population inversion starts 
to break down there can be a large number of longitudinal oscillation modes, 
each producing a separate wavelength in the cavity. Typically, a small number of 
closely spaced wavelengths will predominate, but the resultant beam is not strictly 


monochromatic without special modifications to the laser (see Section 2.5.7). [6] 


2.5.4 Coherence Length 


The preceding discussion leaves the impression that all light generated in the laser 
cavity is coherent. ‘This is not the case. Not all stimulated photons are descended 
from the same spontaneously emitted photon, and the Heisenberg uncertainty 
principle places a lower lint on the precision with which the wavelengths are 
reproduced. So even a laser oscillating in only a single longitudinal mode has 
soine degree of incoherence in its output. 

One measure of the quality of a laser’s output is termed “coherence length’. 
This length is the maximum beam path length within which the beam can coher- 
ently interfere with itself. The practical ramification of this is that the holographic 
system must have all path lengths from laser output to recording planes less than 


the coherence length of the laser. [7] 


2.9.5 Switching 


It. is clear that the amplification properties of the laser cavity are critically depen- 
dent on the alignment and quality of the end mirrors. Misalignment. or surface 
unperfections can destroy the amplification effect and turn the laser into an ex- 
pensive flash bulb. 

This sensitivity can be turned to the laser users advantage. If the mirror prop- 


erties can be accurately controlled then the amplification process can be controlled 
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independently of the pumping process. 

Borrowing a term from electrical engineering, the ability of the cavity to res- 
onate is called it’s ‘Q’. Systems that switch the amplification process on and off 
by controlling the resonance qualities of the cavity are called ‘Q-switched’ or ‘Q- 
spoiled’ lasers. 

Q-switching has a subtle benefit. The population inversion does not occur 
istantaneously at the beginning of the pumping sequence. If the cavity cannot 
oscillate due to low Q, the population inversion can build to a much higher level 
than required for lasing. If the Q is suddenly increased, allowing the cavity to 
oscillate, the peak energy in the pulse can be much higher than if the cavity was 
allowed to oscillate continuously through the pumping sequence. So the power 
available in a Q-switched pulse can be much higher than in an unswitched pulse. |6| 

An additional benefit of Q-switching is the ability to produce multiple pulses 
out of a single pumping process. This allows accurate control of continuously 
pumped lasers, and multiple rapid pulses out of short pulse pumped lasers such 


as ruby lasers. 


2.5.5.1 Electrooptical Q-Switching 


Several physical phenomena allow the electrical modulation of an optical signal. 
The most commonly used in Q-switching lasers is embodied in the Pockels Cell. 

Pockels Cells are formed from a crystal, typically potassium dihydrogen phos- 
phate (KDP). KDP has two distinct optical axes. With an electrical voltage 
applied across the crystal the dex of refraction on one axis is different than the 
other. An incident polarized wave is split into two portions, one portion traveling 
more slowly than the other. ‘The two waves are recombined upon exiting the cell. 
resulting m1 a wave that is elliptically, circularly, or linearly polarized depending 
on the applied voltage. 

The two voltages of interest for Pockels Cell Q-switching are termed the quarter- 
wavelength (A/4) voltage, and the half-wavelength (\/2) voltage. The quarter- 
wavelength voltage will result in a circularly polarized beam at the output of the 


cell. The half-wavelength voltage will result in a wave linearly polarized 90° to 
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Figure 2-5: Pockels Cell Q-switch at (a) quarter-wave and (b) half wave retardation 
voltage. 


the mmput wave. 

The quarter-wavelength retardation system shown in Figure 2-5(a) is a ‘pulse 
off’ Q-switch. During the initial phase of pumping, the A/4 voltage is applied to 
the cell. The output of the Pockels cell is circularly polarized, after reflection the 
circularly polarized beain reenters the Pockels cell and emerges linearly polarized 
90° to the input beam. The polarizer then prevents the beam from reentering 
tle laser rod, spoiling the amplification. When the voltage is removed, the Q is 
restored, and oscillation can start. 

The half-wavelength retardation system shown in Figure 2-5(b) is a ‘pulse on’ 


Q-switch. With no voltage applied the Pockels cell does not alter the polariza- 
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Figure 2-6: Schematic Cavity Dump 


tion of the transmitted beam, and the crossed polarizers prevent oscillation from 
occurring. When switched on the Pockels cell shifts the polarization 90° allowing 
the beam to pass uninterrupted through the two polarizers. 

The half-wave retardation voltage for KDP is 7500V. KDP is also water soluble. 
Normally the crystal is hermetically sealed to prevent exposure to atmospheric 
moisture, but hydrophobic materials and extremely high voltages merit. caution 


when used in a hydrodynamics experiment. 


2.5.5.2 Mechanical Q-Switches 


There are two common methods of mechanically Q-switching a laser. They both 
rely on mechanically rotating a reflective element within the laser cavity. The most 
obvious is to rotate the mirror such that it is only periodically parallel to the fixed 
reflector at the opposite end of the cavity. This system requires extremely critical 
alignment of the shaft upon which the mirror rotates. The alternative to rotating 
a mirror is to use a roof prism, which alleviates the alignment requirement. 
Mechanical switching is not as flexible as electrooptical switching, and is prob- 
ably not practicably capable of implementing the three pulse technique described 


i Section 3.1. 
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Figure 2-7: Cylindrical Transverse Electromagnetic Modes 


2.5.5.3 Cavity Dumping 


Extremely short pulses can be generated using a “cavity dump”. Cavity dumps 
allows the use of two 100% reflectivity end niurrors. When the laser power is 
mlaximized, the cavity dump redirects the beam outside of the cavity, rapidly 
dumping the entire optical energy of the cavity into the output beam. The pulse 
length is reduced to the round trip transit time within the cavity. Cavity dumps 
can be used to switch intermittently and continuously pumped lasers. 

A schematic of a cavity dumped ruby laser is shown in Figure 2-6. Initially. no 
voltage is applied to the Pockels cell. When the flash lamp is fired the horizontally 
polarized wave is passed through the polarizer and not allowed to regenerate. 
Vertically polarized light is kept within the cavity and builds up the vibration mode 
within the cavity. When the energy stored in the Ruby has peaked, the Pockels 
cell is switched to its half-wave retardation voltage. ‘The vertically polarized light 
is rotated to horizontally polarized light by the Pockels cell and passes through the 
polarizer and out of the cavity. The combination of the polarizer, off-axis mirror 
and the Pockels cell essentially forms a voltage controlled mirror, whose reflectivity 


can be change from 100% to 0% and back to 100% during the cycle.[6, 8) 
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2.5.6 Transverse Electromagnetic Modes 


Laser cavities do not only oscillate longitudinally. There are also transverse modes 
of vibration which are very similar to the vibration modes of a drum head. Laser 
specifications typically include a description of which transverse modes are present 
in the output. The Transverse Electromagnetic Modes (TEM) are described using 
the svinbol TEMmng. The first two indices describe the number of transverse 
modes, while the third (often omitted in specification) describes the number of 
longitudinal modes. Figure 2-7 shows cylindrical transverse modes, along with the 
indices describing them. Rectangular modes are also possible. It is not difficult 


to find TEMog ruby lasers suitable for HPIV use. 


2.5.7 Output Frequency Control 


It is possible to dramatically reduce the number of longitudinal modes present in 
the output. The simplest niethod is to use a Fabry-Perot etalon. This is simply 
a very small optical cavity that resonates at precisely the frequency desired in 
the laser cavitv. The etalon effectively spoils the cavity Q for other frequencies, 
preventing them from building up. 

The simplest form of etalon is simply a thin glass plate. Internal reflections 
within the plate selection interfere with the wave to produce a sharp resonance 
peak. The effective thickness of the plate can be changed by rotating the etalon 
about an axis perpendicular to the cavity axis. 

Other etalon designs are possible. Alternate designs include wedge shaped 
plates that can be tuned by translating the plate, as well as designs that form 
an optical cavity using an air gap between two plates. The latter design can 
sometimes be adjusted by varying the gas pressure or temperature within the 


cavity.[6, 8] 


2.5.8 Laser Construction 


Conmnercial lasers for field use are often manufactured in such a way as to limit the 


ability of the user to alter the laser. This increases laser reliability. Commercial 
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lasers for laboratory use are a different matter. As discussed above there are 
many additional components that can be added to the laser cavity to customize 
and control the laser output. Laboratory lasers typically mount. all components, 
including the mirrors and laser rod, on a long rail with additional space for adding 


polarizers, etalons. Q-switches and cavity dumps. 
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Chapter 3 


HPIV System Design 


3.1 


Timing 


The requirements to unambiguously determine speeds less than 10m/s inside a 


volume 5 cm on a side impose severe restrictions: 


Me 


Exposure times must be short enough to ‘stop’ the motion. In a photograph 
the motion allowed during the exposure is related to the feature size of the 
object. During the photographic exposure the smallest features of an object 
should not move more than approximately 1/10th their size or the image will 
be unacceptably blurred. Likewise an acceptable amount of particle motion 
during holographic exposure has been found to be approximately 1/10th of 


a particle diameter|9]. 


The interval between exposures must be short enough to maintain the iden- 
tity of the particles, so tracking and velocity determination is possible. This 
restricts the mean free path length of a particle between exposures to less 
than the distance between particles. There is an additional problein of de- 
termining which particle image is ‘oldest’, so that direction of travel can be 


determined. 


The interval between exposures must be short enough to ensure a majority 
of the particles captured in the first exposure are within the covered volume 


during the following exposures. 
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For a 10jm particle, these requirements become: 
1. A ls exposure results in 1/10th diameter of travel during the exposure. 


2. A 10us exposure interval results in 10 particle dianieters of travel between 


exposures at 10in/s. 


The exposure interval also places an upper limit on the concentration of seed 
particles. Assuring the particle mean free path between exposures is less than the 


average distance between particles requires that 


particles particles 


| 5 6 
jo T(100jan)* = 2x NO ce 4x 10 3 
Where p, is the seed particle density, and the mean free path of the particles is 
ten particle diameters, 1O0Qum. 

The direction problem is subtler. Keeping the average inter-particle distance 
greater than the mean free path length between exposures significantly simplifies 
the task of identifying pairs of particle images as being the same particle. Unfortu- 
nately, even when the two images of a given particle have been identified, there 1s a 
180° ainbiguity in the velocity of the particle. A simple solution to this problem is 
to use three pulses rather than two, with differing time intervals between pairs of 
pulses. For exainple, the second ls pulse is sent Sys after the first, then the third 
is sent 10us after the second. The pulse intervals are sufficiently short to permit 
particle identification, while the difference in intervals allows determination of the 
tine order of the images. An additional advantage of the triple exposure tech- 
nique could be estimation of fluid acceleration. This technique assumes the flow 
speed will not change appreciably over the distances traveled by particle during 
the exposure sequence. If the particle were to radically decelerate to less than half 
of the speed it traveled during the first set of exposures, it would be possible to 


confuse the tine order of the exposures. 
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3.2 Recording and Digitization 


3.2.1 CCD Requirements 


The ultimate goal for computationally reconstructed HPIV is a completely digital 
process for determining the flow field throughout the volume. The most likely can- 
didate for the ‘digital recording plane’ is a Charge Coupled Device (CCD) Array. 
Unfortunately. the resolution of CCD arrays is significantly less than holographic 


filnin. 


3.2.1.1 Characteristics of CCD arrays 


A CCD Atrravy is a semiconductor architecture that transfers charge through stor- 
age areas. The CCD array is essentially an array of charge “buckets”. Photons 
are captured by the buckets and converted to charge, the buckets only perform 
this conversion during a controllable integration period. More photons, or higher 
energy photons, hitting the array during the integration period result in a higher 
element charge. After integration the CCD Array support electronics scan through 
each bucket and report the charge level in each bucket. This charge is proportional 
to the intensity of the image focused on the array./10) 

For holographic applications the fundamental characteristic of CCD Arrays 
that differentiates them from film is the limited number of recording elements. 
High quality photographic film can have grains smaller than lym over an area 12- 
20 centimeters square, or better than 120.000 120,000 recording elements. The 
best CCD arrays currently commercially available have element sizes on the order 
of 10jum, but the total array size is only about 5 centimeters on a side, yielding 
all array of approximately 5000x5000 elements. This is the root. of the limitation 


in holographic applications. 


3.2.1.2 CCD Resolution Requirements 


3.2.1.2.1 Diffraction patterns of circular objects Parrent and Thomp- 


son analytically derived the diffraction pattern for an opaque circular object in a 
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coherent field. [11] For a circular particle, the field and intensity, respectively. are: 


ka? - kr? kar 
Wir) = coskz — —sin[kz+ oa Ay 
1 = ae 


Aw 


ka? kr? Or 
+i sinks + “oo [e+ 5 Ja ( )| 
in ba z 


2ka? — kr? ka ka’ kar \]° 
he |e ee sin SA ( =| — A (“")] (3.1) 
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Where, r is measured from the center of the object projected into the recording 
plane, and A(a@) = J;(a@)/a, a Bessel function of the first kind divided by its 
argument. 

An illustrative example is plotted in Figures 3-] and 3-2. The parameters 
more characteristic of HPIV systems result in plots that do not acceptably show 


the salient features of the mtensity distribution. 


3.2.1.2.2 Recording the Diffraction Pattern Equation (3.1) describes the 
diffraction pattern recorded for a single particle with circular cross section. All 
previous discussion has assumed this diffraction pattern is perfectly recorded. In 
reality the recording medium has limited bandwidth (7. e. a grain), and limited 
physical extent. 

Figure 3-1 shows a high frequency sine function modulated by a lower frequency 
A(a) envelope function. The envelope function is dependent on both particle 
diameter and particle distance from the recording plane. The high frequency 
sine function is dependent only on particle distance. Accurately recording the 
high frequency component is the most critical factor in obtaining good depth 
discrimination. 

Vikram|3] recommends recording the first three side lobes formed by the bessel 


envelope function. The spacing of the fine fringes, Ax, in the third lobe can be 
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Figure 3-2: Enlarged plot of the third side lobe in Figure 3-1 


found using the argument of the sine term in (3.1), kr?/2z: 


hae Gis 5 
eI a Paes = ik 
22 ae 
Te We az 
r 
i es (3.2) 
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Equation (3.2) should be evaluated in the region containing the third side lobe. 


The mid-radius of the third side lobe, 7s is: 


OT) 07 = Cl ye? oe 


Qn) Coc U0, LOM edtoo7s, 1o4t70. 





kary 
ook 
= = r, + Ci) _ z | a 
i 9 aa ey (ay, is a) 
rf; =02 (es fos) 17 eo (3.3) 
: oka * ; ha . 


Equations (3.3) and (3.2) give 7; = 0.189cm and Ar = 3.65 x 107%cm for 
the parameters plotted in Figures 3-1 and 3-2. Usimg z = 5cm and a = 10pm. 
with A = 690nm, more reasonable for HPIV purposes, 73 = 0.472cm and Ar = 
ol x 10F emt, 

The parameter 73 dictates the minimuni size of the hologram, with no mag- 


nification, and is easily met. The parameter Az dictates the minimum grain of 
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the hologram recording medium. Since enlargement or reduction of the hologram 
during recording is easily accomplished using lenses, the absolute size of the fringe 
spacing is not important. Rather, the important quantity is the ratio the fringe 
spacing to the total size of the system, so that the fine fringes produced by par- 
ticles at opposite extremes of the volume can be resolved. The sampling theorem 
suggests that the recording medium should be able to sample at twice the high- 
est. spatial frequency required, or 1/2 the fringe spacing. To achieve the desired 


accuracy. a CCD array would have to have 


2-0.0Sscem 


——————_——. = 13680 pixels square. 
7.31 x 10-4cm epee 


This resolution will adequately record the particle diffraction patterns. Unfortu- 
nately no CCD arrays are commercially available that reach this resolution. The 


best currently available are in the range of 5000 pixels square. 


3.2.2 Exposure Requirements 


While the performance characteristics of CCD arrays are well known it is difficult 
to predict the amount of power incident on the array from the suspended particles. 

The resolution analysis was performed assuming the finest necessary fringes 
were detected. The incident power depends heavily on the power output of the 
laser, the absorption of the fluid, and the reflectivity of the particles, and the 
viewing window material. 

The exposure problem is different in character from: the resolution problem. 
The resolution problem has to wait for better CCD technology, or lower the ac- 
curacy requirements for HPIV significantly. The exposure problem can be solved 
many ways, the simplest (but not necessarily the least expensive) way is to increase 


the laser power. 
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Chapter 4 


Simulation 


4.1 Motivation 


A primary factor preventing HPIV from becoming an easily and widely used tech- 
nique is the data reduction effort. As previously mentioned, physical reconstruc- 
tion techniques are simple, but require costly equipment, and inordinate amounts 
of time. If coniputational reconstruction methods can be perfected, HPIV could 
be become more practical. 

There are several possibilities for computational reconstruction algorithins. 
one of which (Onural’s method) is presented in Section 5.1. The methods ust 
be evaluated with realistic input in order to gauge their effectiveness. Obtaining 
realistic input 1s difficult. There have been experiments done using holographic 
techinques to study small ocean biological specimens, but obtaining the actual 
holograms produced is difficult. An additional difficulty is that you don’t know 
precisely what was contained in the original image so calibrating the reconstruction 


algorithm is nearly impossible. 


4.2 Hologram Synthesis 


There are several methods applicable to synthesizing holograms similar to those 
expected from HPIV. The most straightforward methods calculate the individual 


contribution from each particle to each pixel on the simulated recording plane. 
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These methods are the most accurate and flexible but suffer from excessive com- 
putational requirements when large numbers of particles are being simulated. 


One way around this problem is hinted at by Equation (2.18): 


o(a,y) =) Un(2,y) @ hi, (2,y) 


n 


We can subdivide the regiou being simulated into a large number of thin slabs. 
Within each slab there can be a large number of particles, but each will be treated 
as if it had the saine z coordinate. We then calculate the contribution of each 
slab to the total hologram rather than the contribution from each particle. With 
a large number of particles we cau achieve significant, performance increases with 
little degradation of quality. 


Evaluating (2.18) is quite efficient. Using the convolution theorem: 


bey) = > Un(a,y) @h,, @,y) 


> FoF (Un (x,y) F [h:, (2, y)] 
= Fl So Un Ree al (Gee cs) (4.1) 


“n 


Where. (a2,W,)— 7 |U (2,9) and 712, (u,,0,)— 7 [hn (ay), But Hoes, 


can be analytically evaluated [12], 





es 9 9 
FiO) eo |i a| sex) i ie (w*. ++ %)| (4.2) 
yielding. 
U (0.9) = FOS Uy (Wes ey) exp [j= (w? + 02) (1.3) 
‘ = Ar Ir y 


The algorithm to evaluate (4.3) is shown schematically in 4-1, and computer 


code to perform the algorithm is given in Sections A.1.1 and A.1.2. 
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Figure 4-1: Hologram synthesis algorithm 
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4.3. HoloSynth Simulator 
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Figure 4-2: HoloSynth user interface 


4.3.1 Introduction 


HoloSynth assists in synthesizing holograms similar to those obtained from Holo- 
graphic Particle Image Velocimetry. It provides the user with a simple way to 
place any number of particles into a cubic volume in preparation for generating 
a synthetic holograin. The user can manually place particles anywhere in the 3- 
dimensional region as well as fill arbitrary rectangular solid regions with randoinly 
distributed particles. 


HoloSynth uses the method for synthesizing holograms described in Section 
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Figure 4-3: HoloSynth main controls 


4.2. In addition it prepares an mput file for an alternate hologram synthesizer 
so that different. synthesis algorithms may be easily compared, without modifving 
the HoloSynth source code. The output file is also useful in analyzing the results 
of reconstruction algorithms. 


HoloSynth was developed under Red Hat(Q)Linux version 5.1 and tested under 


Microsoft Windows NT@®). 


4.3.2 The Interface 
4.3.2.1 The Main Controls 


Figure 4-3 shows the main controls for HoloSynth. 
Each button will open up a new window that controls a particular aspect of 


HoloSynth. The windows are: 


1. 3-D View of the current particle configuration. 


bo 


Plane View, which allows manual placement of particles on an arbitrary 


plane. 


3. Select Region, which allows selection of an arbitrary rectangular solid region 


withm the view space. From this control panel you can: 


e Clear all or any part of the volume. 
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o Cesar 1.0: 3-D View 





Figure 4-4: 3-D View window 


e Fill all or any part of the volume with randomly placed particles. 


4. The Synthesizer Settings window which controls the physical characteristics 


of the simulation used for the hologram synthesis. 


Through the file menu the user can open and save particle fields, and exit 


Holosynth. 


4.3.2.2 The 3-D View 


The 3-D View. shown in Figure 4-4, provides a three dimensional view of the 
current particle configuration. ‘There are three basic functions controlled by this 
window: the magnification (zoom) of the view, the angle and translation of the 
view, and the position of the active and two clipping planes. No direct editing 
takes place in this window, it is merely a visualization aid. 

The viewing “camera position” is controlled by the rollers and shders to the left 
and below the viewing window. The slider above the viewing window controlling 
the magnification of the camera. 

‘The various view control planes are controlled using the slide bars to the right of 


the 3-D view window. They appear as colored squares drawn around the periphery 
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Figure 4-5: HoloSynth plane editor 


of the cube. Only the active plane is visible in Figure 4-4, the clipping planes are 
hidden at the ends of their travel. On a color display the active plane is bright 
ereen. the rear cluping plane is blue, and the front clipping plane is magenta. The 
active plane determines the z-coordinate that the plane editor uses when adding 
lew particles. The front and rear clipping plane restrict the range of z coordinates 
visible in both the plane view and the 3-D view. In all cases, particles lying on 
the active plane are visible within both views. If the user wants to view only the 
particles on the active plane, then he places both clipping planes at the same limit 


of travel. 


4.3.2.3 The Plane View 


The plane view projects all particles visible in the 3-D view onto a single plane. 
Particles that are on the active (green) plane in the 3-D view are highlighted green 
in the plane view. Particles can be manually positioned or deleted on the active 
plane. 

‘To add a particle, simply left click within the black portion of the active plane 


view. Accurate positioning of particles can be performed in two ways: 
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o HoloSynth 1.0: 3-D Selection 
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Figure 4-6: 3-D region selection 


1. Turn on the coordinate grid. At the bottom of the plane view window is a 
selection list allowing presentation of grids of various densitites. High grid 
densities are convenient at high magnifications. While a grid is displayed 
the particles automatically snap to the nearest grid intersection. This can 
be controlled by the Display->Snap to Grid menu item, and can also be 


toggled using the ‘g’ key. 


No 


Use the real-time coordinate diplay. Any time the mouse is over the edit 


window its complete 3-D coordinate is show at the bottom of the window. 


Deletion is accomplished by dragging with the right mouse button over the 
region containing the particles to be deleted. A red selection border will appear 
showing the selected region. ‘To delete inside the selected region, use the menu 
item, Edit->Delete Inside Selection or the delete key. Using the Edit->Delete 
Outside Selection menu item will clear all particles visible in the plane editor 
but outside the selection box. The region deleted is bounded in the x and y di- 
rections by the selection box, and in the z dimension by the clipping planes. If 
a particle is visible in the plane editor, regardless of whether it is on the active 


plane, it can be deleted. 
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4.3.2.4 Select Regions 


The 3-D Selection window, shown in Figure 4-6, allows the user to define a solid 


rectangular sub-volume within the 3-D view and perform the following operations: 


1. Clear the entire volume. 


NS) 


_ Clear all particles inside the selected sub-volume 
3. Clear all particles outside the selected sub-volume 


4. Fill the entire volume with randomly placed particles. The slider at the right 


of the window controls how inany particles will be added. 


5. Fill the selected sub-vohume with randomly placed particles. The shder at 


the right of the window controls how many particles will be added. 


6. Fill outside the selected sub-volume with randomly placed particles. The 


slider at the right of the window controls how many particles will be added. 


The sub-volume is displayed as a translucent yellow box with the 3-D view. Each 
wall of the selection region is controlled by its own slider. The region will only be 


displayed in the 3-D view if the Show Selection Box button has been selected. 


4.3.2.5 Hologram Synthesis 


Figure 4-7 shows the synthesizer controls. This allows easy control of the physical 
size of the simulation as well as recording plane clistance and number of planes used 
to subdivide the simulation volume. The generated hologram will be displayed in a 
hew window when the computations have completed. If the product is satisfactory, 
the user may save the hologram to a file using the Truevision Targa format or the 


Adobe Photshop Raw format. 


4.3.3 Requirements 


HoloSynth is written completely in C++. The Graphical User Interface uses 


the freely available platform mdependent user interface toolkit FLTK{13]. The 
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Figure 4-7: Synthesizer settings 
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hologram synthesizer use the FF TW fast fourier transform hbrary available freely 
from MIT/14]. The program runs under X-Windows, and Microsoft Windows with 


OpenGL/15j or the Mesa 3-D graphics library (an OpenGL functional equivalent ){16]. 


4.3.4 Program Architecture 


At its heart HoloSynth is a fancy list editor, editing a list of particle coordinates 
and radu. The core data structure is, appropriately, a linked list of particle objects. 
This allows an arbitrary number of particles to be added, restricted only by the 
inemory available. 

Figure 4-8 shows the schematic relationship between all program modules that 
make up HoloSynth. The complete source code to HoloSynth is presented in 


Appendix A.1. 
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Figure 4-8: HoloSynth architecture 
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Chapter 5 


Computational Reconstruction 


5.1 Onural’s Method 


In 1985, Levent. Onural proposed a technique based on a truncated inverse filter|12, 
17}. The method is summarized here, source code is presented to perform the 
algorithm in Sections A.2.1 and A.2.2. 


Starting from (2.17), 
v2(t,y) = U(x, y) @ h.(a, y) 


We let U (2, y) = 1—a(a,y) where a(z,y) is the opacity of the given plane. Now 


write the recording plane intensity, 


I(z,y) = v2 (2,y) v3 (2,y) 
= |(1—a(zx,y)) @h, (2, y)|° 
= l—-a'‘(z,y) @hj (z,y) —a(z,y) Oh, (z,y) 
+ a(x, y) @hz(x.y)}" (5.1) 


= l—a(z,y) @ th, (a, y) +h, (2, y)] 


qn 
ho 
SS 


= 1-—a(z,y) ® 2Re {h,; (2, y)} (5. 


Where |a (2, y) @ hz (x, y)|° was dropped as insignificant in (5.1). Equation i.) 
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shows the intensity modeled as a DC shifted linear system with mpulse response. 


g(x,y) = 2Re th. (x, y)} (5.3) 


No inverse of g(z,y) in the normal sense exists, However, Onural showed in 
[12] that a series approximation to the Fourier transform of g can be found. If 
Cio =F olay || hein, 

ie est 


: ee Mf — 1 — [log, k] _ 
Gi} (eran) eon ge (Has) + fe LN (wren) (04) 


where {log, A] is the integer obtained by truncating the value of log, & and, 


Az . 
Sone. (ne cos = (2k + 1) (we + ws) (9.5) 
yielding. 
Cecds CA et aed reat ic) (5.6) 


Looking back at the holograin recording given by (5.2), 


Ota) CORE ay) = a i) 
F la(z,y)|G (wr wy) = F[l—I(z,y) 
an(z,y) = F  |F[l—1 Gy) Gwe w2.4,)| (6.7) 


vielding the original opacity function. 


5.2 Simulation Results 


5.2.1 Synthesis and Computational Reconstruction 


Figure 5-2 is a portion of a hologram generated by HoloSynth. The particle config- 


uration consisted of 2328 particles. 2000 of the particles were randomly distributed 
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Figure 5-1: Schematic of the particle configuration used to generate Figure 5-2 


throughout a 5em cube. Along the center plane parallel to the wavefront 328 par- 
ticles were positioned to spell the letters “MIT”. Figure 5-1 illustrates a schematic 
of that configuration. The synthetic hologram itself was 2048 x 2048 pixels, and 
the volume was subdivided into 1024 planes, for a z-axis discretization of 0.04mm. 
The simulated particle diameter was 25um. The portion of the hologram shown m 
5-2 is 512512 pixels and its upper left corner is 205 pixels from the left extremity 
of the origmal and 741 pixels from the top. Although it contains the ‘M’ from 
“MIT” , it is impossible to distinguish from the diffraction patterns produced by 
the other particles. 

Figure 5-3 through 5-8 show the results of reconstructing Figure 5-2 at various 
distances between 30.0cm and 30.5cm from the recording plane using Onural’s 
mnethod with AJ = 5. The images generated are the opacity functions, with black 
bemg transparent, white being opaque. 

In Figure 5-3(a) the ‘M’ is clearly visible on a computer monitor with some 
slight clouds from out-of-focus particles. Figure 5-4 is the same hologram recon- 
structed at 30.1cm. Here there is little visible difference from the previous figure. 


This is an artifact of the way the decoder transforin real valued intensities into 
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gray scale bitmap images. 

To translate from an arbitrary field of real number to a field of 8 bit integers. 
the decoder first finds the maximum and minimum values present in the real array. 
it then scales and translates each of these values to fill the range from O to 256. 
in order to maxinize use of dynamic range available in an 8 bit gray scale image. 

The 30.lem reconstruction has a much lower maximum intensity than the 
30.0 cin reconstruction. Volume reconstruction, where many planes will be recon- 
structed at a time, in order to reconstruct flow fields, will have to keep track of 
this scahng issue between planes. 

The important thing to note about this series of decodes, is that while the 
source hologram was constructed using a z discretization of 0.04mm, the recon- 
struction algorithm z-discrimination was of the order of millimeters. Thus, the 
original image could have been discretized into 64 steps, or 0.78mm slabs, requiring 


1/16th the computation time with little or no loss of reconstruction accuracy. 


5.2.2 Image Postprocessing 


The next step in reconstruction is to automatically detect particles in each plane 
and register their z, y, and z coordinates, rather than simply displaying an im- 
age. Repeatedly reconstructing at different depths and detecting the particles will 
replace the time consuming 3-D scanning. 

In Figures 5-3 through 5-8 subfigure (b) and (c) show the results of applying 
some simple image processing techniques to the original reconstructions to enhance 
the output. Subfigure (b) in each figure applies ‘edge detection’, which simply 
filters the image marking any region where the gradient of the image exceeds 
a certain threshhold. In the case of the grayscale images being used here, the 
gradient. is sumply the difference in grayscale value between adjacent. pixels.{18] 
Subfigure (c) in each figure ‘threshholds’ the image, changing a pixel to white 
(value 255) or black (value 0) depending on whether or not the pixel value exceeds 


a preset. threshhold.? 


'The image processing was accomplished with the freely available GNU Image Manipulation 
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The combination of edge detection and threshholding can be used very effec- 
tively to automatically detect particle positions. At the end of the image process- 
ing, there are only regions of black (value 0) with very small regions of white (value 
255). The white regions correspond to particles and can easily be automatically 


tabulated and mapped into x and y coordinates. 


Program (GIMP, http://www.gimp.org), which is very similar to Adobe Photoshop®. 
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Figure 5-2: Portion of a synthetic hologram. The original synthetic hologram was 
2048 x 2048 pixels. ‘This portion is 512512 pixels, cropped from the orig- 
inal hologram to make the images presentable in printed form. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de- 
tection. 





(c) Image after enhancing with edge de- 
tection and threshholding. 


Figure 5-3: Figure 5-2 decoded at a depth of 0.300m. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de- 
tection. 





(c) Image after enhancing with edge de- 
tection and threshholding. 


Figure 5-4: Figure 5-2 decoded at a depth of 0.301m. 
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(b) Image after enhancing with edge de- 
tection. 


(a) Raw Reconstructed Image 





(c) Image after enhancing with edge de- 
tection and threshholding. 


Figure 5-5: Figure 5-2 decoded at a depth of 0.302m. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de- 
tection. 





(c) Image after enhancing with edge de- 
tection and threshholding. 


Figure 5-6: Figure 5-2 decoded at. a depth of 0.303m. 
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(b) Image after enhancing with edge de- 
tection. 


(a) Raw Reconstructed Image 





(c) Image after enhancing with edge de- 
tection and threshholding. 


Figure 5-7: Figure 5-2 decoded at a depth of 0.304m 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de- 
tection. 





(c) Image after enhancing with edge de- 
tection and threshholding. 


Figure 5-8: Figure 5-2 decoded at a depth of 0.305m. 
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Chapter 6 


Conclusions and 


Recommendations 


It is clear that computational reconstruction of HPIV is practicable. ‘The proce- 
dure that Onural proposed in 1985 was not then practical due to the computational 
power cominonly available. Modern computing can easily handle realistic holo- 
grams and the numeric techniques of normal PIV could be extended to deal with 
the particle identification problem in HPIV. The simulator presented in this thesis 
will be useful in any setting when attempting to validate particular reconstruction 
algorithms. 

Future work should concentrate on extending the reconstruction prosses fully 
into three dimensions and developing particle identification techniques. Automat- 


ically determining the velocity vector for a particle is the ultimate goal. 
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Appendix A 


Source code 


A.1 HoloSynth: Hologram Synthesizer. 


This is the source code for HoloSynth, a program for easily creating particle configuration 


and holograms for use calibrating computational reconstruction algorithms. 


A.1.1 OnuralSynthesizer.hpp 


#ifndef ONURALSYNTHESIZER_H 
//linclude guard 
#define ONURALSYNTHESIZER-H 1 


#include “Synthesizer .hpp” 
#include "ParticleList.hpp’ 
#tifndef FLMATH 

#define FLMATH 

#include <FL/math.h> 
#endif 

#include "RGBImage.hpp’” 


/*\ class OnuralSynthesizer 
* 
*\ brief The algorithm use for actually generating the hologram 
* 
«The forward diffraction algorithm implement a method described in 
* “Digital Decoding of In-Line Holograms” a PhD thesis by Leven Onural 


7 


class OnuralSynthesizer:public Synthesizer 


{ 
public: 
//! created with a reference to the current particle list 
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OnuralSynthesizer(ParticleList *pl); 

~OnuralSynthesizer(); 

//! Generate, organized all the FFTs and diffraction calculations. 
int generate(): 


private: 
RGBlnage *lologram: 


void diffract(fftw_complex *image); 

void clearField(fftw_complex *field, fftw_real val); 
void plotParticle(particle «c.fftw_-complex *plane); 
void plot(int p_X, int p_Y, fftw-complex plane): 
void addto(fftw_complex *al, fftw_complex *a2); 
void MagtoRGB(fftw_complex *data); 


void cleanup(): 


fftwnd_plan fftwFwdPlan:; 
fftwnd_plan fftwBwdPlan:; 


fftw_complex *oplane: 
fftw_complex *hplane; 


#endif 


A.1.2 OnuralSynthesizer.hpp 


#include <FL/F1.H> 

#include ’ProgressUI .hpp” 
#include “ScrollImageUI .hpp> 
#include <FL/F1l_Image.H> 
#include "OnuralSynthesizer.hpp” 


#include "StatusBox.hpp” 
#define MIN(a, b) (((a) < (b)) ? (a) : (b) 
#define MAX(a. b) (((a) > (b)) ? (a): (b)) 


OnuralSynthesizer::OnuralSynthesizer(ParticleList. «pl) 
‘Synthesizer(pl) 


} 


OnuralSynthesizer::~OnuralSynthesizer() 


{ 


cleanup(): 
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} 


int OnuralSyuthesizer::generate() 


{ 


int progresscount=0; 


progressUI—>set Max(5+5*d_Planes); 
progressU]>update(Q); 
progressUI— reset (); 
progressUI—show(); 

Fl::wait(0.0); 


particle *curr; 
int pixels, 11, pcount, Ncalc, Mcalc; 


fftw_real zgrain, Zmin, Zmax; 


OOOO ORO ORO ORI ORI OOO IORI OCI IORI IIR IARI IK ICE 


* Initialization. 


FOO OOOO OOOO IS IORI ISI II II IC Ka / 
pixels=2*(2*N /2+1)*2*M; 


Ncalc=2*N; //Double the size for padding out the edge effects. 
Nealc=2eM- 


if(progressU]—setStatus("Allocating memory...")){ 
cleanup(); 
progressUI—hide(); 
return 0; 


}: 


hplane=(fftw_complex *) calloc( pixels,sizeof(fftw_complex)); 
oplane=(fftw_complex *) malloc( pixels*xsizeof(fftw_complex)); 


progressUI-update(++progresscount ); 

zerain=(fftw_real)1/(2*d_Planes);//The thickness of a HALF plane 

if(progressUIl—setStatus("Generating FFTW Plans...")){ 
cleanup(): 


progressUI—hide(); 
return 0; 


fftwFwdPlan=fftw2d_create_plan(Ncalc, Mcalc, FFTW_FORWARD. 
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FFTW_IN_PLACE); 
fftwBwdPlan=fitw2d_create_plan(Ncalc. Mcalc, FFTW_-BACK WARD, 
FFTW_IN_PLACE); 


progressUI > update(++progresscount ); 
//Remember we will eventually need to split out the vector end of things. 


pList—sortListOnZ(); 


curr=pList—first(): 


for (il=1;i11<=d_Planes;i1++){ 


Fl::wait(0): 


if(progressUI—setStatus("Clearing Object Plane...")){ 
cleanup():; 
progressUI— hide(); 
return 0: 


: 


clearField(oplane,1.0); 
progressUI— update(++progresscount); 


Zmin=(il-1)*2*zgrain-0.5:; 
Z=(Zmin+zegrain)*depth;// Physical Z 
Zmax=Zmin+2*zgrain; 


Deount—0:- 


if(progressUI—setStatus("Plotting Particles...")){ 
cleanup(); 
progressUI—hide(); 
return 0: 


ie 


while((curr!=NULL)&&(curr>Z>Zmin)&&(curr>Z<=Zmax)){ 
plotParticle(curr.oplane): 
peount+-4; 
CURT —CUrr— ext: 

} 


progressUI— update(++progresscount): 


if(pcount>0){ 


if( progressUI—setStatus("Performing Forward FFT...")){ 
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cleanup(); 
progressUI>hide(); 
return 0): 


}; 


fftwnd_one(fftwFwdPlan, oplane, 0); 
progress UI update(++ progresscount); 


if(progressUI—setStatus("Diffracting...")){ 
cleanup(); 
progressUI—hide(); 
return (0); 
ts 
diffract(oplane); 
progressUI>update(++progresscount ); 


if(progressUIsctStatus("Adding...")){ 
cleanup(); 
progressUI—hide(); 
return 0; 
i 
addto(hplane,oplane): 
progressUI— update(++progresscount); 
} 
progresscount=9%*11+ 2; 
progressUI— update(progresscount ); 


} 


hplane[0].re= 1000000; 


progressUI—setStatus("Performing Backward FFT..."); 
progressUI—> update(++progresscount ); 
fftwnd_one(fftwBwdPlan, hplane, 0); 
progressUI—setStatus("Calculating Intensity..."); 


MagtoRGB(hplane); 

//The hologram is finished, now display it... 

Fl_Image *image=new F1_Image(hologram—data(), N, M); 
progressUI—setStatus("Displaying..."); 
progressUI—update(++progresscount); 


ScrolllmageUI *si=new ScrolllmageUlI (image , 0,0, 500,500); 


si—show(); 


progressUI—hide(): 
cleanup(); 


return |; 
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void OnuralSynthesizer: 


{ 


-addto(fftw_complex *a, fftw complex *b) 
int il; 
for (i1=0; il < (4 * N*M ); i1++){ 
if( (fabs(b[il].re) >=0.0000001)& & (fabs(b[i1].i1m) >=0.0000001)) 


alil|.re += blil}-re; 
alil].in. += bfil].im; 


} 


void OnuralSynthesizer::clearField(fftw_complex *field, fftw_real val) 


{ 


int il; 
for(il=0; il < (4* N*M); i1++){ 
field[il].re = val; 
field{il].im = 0.0; 


} 
} 
void OnuralSynthesizer::plotParticle(particle *c,fftw_complex *plane) 
{ 


int x, y, 1; 


x=(int)((c>X + 0.5)«N); 
y=(int)((0.5 - c+ Y)*M): 
r=(int)(c3RxN); 


switch (r) 


{ 


case 6: 


plot(x-1, y+1, plane): 
plot(x-1, y-1, plane); 

plot(x+1, y-1, plane); 
plot(x+1, y+1, plane); 
case 2: 


case 3: 
case 4: 
case 5: 


plot(x, y-1, plane); 
plot(x, y+1, plane): 
plot(x-1, y, plane); 
plot(x+1, y, plane); 
default: 


plot(x, y, plane); 
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} 


void OnuralSynthesizer::plot(int p-X, int p_Y, fftw_complex »*plane) 
{ 
iIntapos=23° p.1 + IN pak: 
if(pos>0){ 
plane[ pos].re=0.0; 
plane! pos].im=0.0; 


} 


#define F_PI 4.0*3.1415926535 
void OnuralSynthesizer::diffract (fftw_complex *image) 
{ 

int x, y, Ncalc, Mcale; 

long ik: 

fitwereal 1Z.-ce: 

fftw_real wx, wy, a, Ca, Sa; 


INGale=2+1N- 
Meale=2+hi- 


lz=L+*(distance-Z); 
ike): 


/x FFTW puts frequency space origin at the upper right corner of the 
* array. Positive frequencies build from the upper left corner 
* diagonally down through the upper left quadrant. The other corners 
* of the array the array are ALSO at the origin, so the very center 
* of the array 1s the HIGHEST frequency. 


¥* 


In order to use the diffraction filters we MUST take that into 
account, that is what the strange trigraph statements take care of 


¥* 


¥* 


below */ 
for (y=0;y<Mcalc;y++) 


{ 
wy = (y<M) ? y/height : (y-Mcalc) /height; 
wy = (y<M) ? y/height : (Mcalc-y)/height; 
wy * = wy: 
for(x=0;x< Neale;x++){ 


1) 


ne 


wx=(x<N) ? x/width : (x-Nealc)/width; 
// wr=(r<N) ? /width : (Nealc-x)/width; 

a=lz*(wy+wx*wx)*3.1415926535; 
a=l2x (wyt+ weewe)/PF_ PI; 


— 


i, 
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Ca=COs( a); 

sa=sin(a); 

cc=imagelik].re*ca-imagelik].im-*sa: 
image[ik].im=image[ik].rexsa+image[ik].im*ca; 
image[ik].re=cc; 

ik++; 


} 


void OnuralSynthesizer::MagtoRGB(fftw_complex *data) 
int i, j. ij, grey, Nealc, Mcalc: 
fftw_real v, vmax, vmin, vl. v2: 


Nealc=Z24N:- 
Mcalc=2*M: 


hologran=new RGBImage(N,M); 


vinax=0.0; 
vmin=le6: 
=U: 


//Scale the result to take maximum advantage of 8bit range 
for(i=0:1< Mcalce*Nealce;it++){ 

v1l=datali].re; 

v2=datali].im: 

v=sqrt(vl*vl+v2«v2); 


vmax=MAX(v, vmax); 
vmin=MIN(v, vinin); 


} 


vmax=1/(vmax-vmin); 


RGE Color c: 
for (i=0;i<M;i++){ 
for(j=0;j<N;j++){ 
1j=1* Ncalce+]; 


vl=datalij].im; 

v2=datalij].re; 
v=(sqrt(vl*vl+v2*v2)-vmin)*vmax; 
datalij].re=v; 





A.l. 


} 


grey=(int )floor((fftw_real)(v*255.0) ): 
c. B=grey; 
c.R=grey; 
c.G=grey; 
hologram—plot(j, 1, ¢); 


void OnuralSynthesizer::cleanup() 


{ 


free(hplane); 
free(oplane); 
fftwnd_destroy_plan(fftwFwdPlan); 
fftwnd_destroy_plan(fftwBwdPlan); 
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A.2 Decoder: Onural Reconstruction Algorithm. 
This is the source code for an Onural reconstruction algorithm. 


A.2.1 OnuralDecoder.hpp 


#tifndef ONURALDECODER-_HPP 
#define ONURALDECODER-HPP 1 
#include "Decoder .hpp’ 
#include ’ProgressUI .hpp™ 
#Hinclude <fftw.h> 
class OnuralDecoder: public Decoder 
{ 
public: 
OnuralDecoder(); 
virtual uchar* decode(uchar *data, int N, int M, 
float W, float H, float L, 
float D); 


private: 
int d_ImageN; 
int d_CalcN; 
int d_ImageM; 
int d_CalcM:; 


float d_tImageW; 
float d_ImageH; 

float d_ImageLambda: 
float d_decodeDist: 


fftw_complex *d_Field; 


Wears) coal: 


void OnuralDecoder::noDClquad(fftw_complex *field, fftw_real average Value); 
fftw_real OnuralDecoder::aprinv(fftw_real zz, int Mstage); 
void OnuralDecoder::onural(fftw_coniplex *hg, 

int N, int M, 

fftw_real L, fftw_real Z, 

fftw_real X, fltw_real Y.int Mstage); 


fftw_real Magnitude(fftw_complex »*field); 


void uchar2field(); 
uchar *feld2uchar(): 
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fftw_real *clipField(): 


uchar *OnuralDecoder::mainLoop(); 
ProgressUI *progress; 


#endif 


A.2.2 OnuralDecoder.cpp 


#include "OnuralDecoder.hpp’ 
#include <fftw.h> 

#include "StatusBox.hpp” 
#include <FL/F1.H> 

#include <math.h> 


#define PI 3.14159265359 

#tdefine TP! 6.28318530718 

#define MIN(a, b) (((a) < (b)) ? (a): (b)) 
#define MAX(a. b) (((a) > (b)) ? (a): (b)) 


OnuralDecoder::OnuralDecoder() 
{ 

d_ImageN=0; 

d_ImageM=0; 


d_ImageW=0: 
d_lhmageH=0; 
d_InnageLambda=0; 
d_decodeDist=0: 
eedata—0; 

progress=new ProgressUI(); 


j 


TO 


uchar *OnuralDecoder::decode(uchar* data, int N. int M, float W. float H, float L. 


float D) 


d_ImageN=N; 
deCaleN=22N- 


d_ImageM=M: 
d.CaleM=2+M, 


d_ImageW=W; 
d_ImageH=H; 
dtImageLambda=L: 
d_decodeDist=D: 
d_data=data: 
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return mainLoop(): 


/* Remove the DC component from only the first Quadrant, the other 
* three quadrants are ALREADY ZEROx / 


void OnuralDecoder::noDClquad(fftw_complex *field, fftw_real average Value) 
te 
int 1, J, pos: 
for(i=0; i<d_ImageM; i++){ 
for(j=0; j<d_ImageN; j++){ 
pos=ixd_CalcN-4): 
field[pos].re -= average Value: 


} 


fftw_real OnuralDecoder::aprinv(fftw_real zz, int Mstage) 
fftw_real cc,ss,powKa; 
int k,powk: 


eC =c0s|((27); 


Ssa=l): 


if( Mstage>1){ 
SS=— xs 272): 
cc+=ssxsin(2*zz): 
ss* =2*cos(2*zz); 

} 

for (k=3;k<=Mstage;k++){ 
powK=4<(k-3); // powK=2«x(k-1), fine for small k, and FAST 


powKa=(fftw_real) powK «zz: 


cc+=ss*sin(powKa); 
ss* =2*cos(powKa); 


} 


return 2xcc/Mstage; 


j 


void OnuralDecoder::onural(fftw complex *hg, 
int N, int M., 
ffitw_real L, fitw_real Z, 
fftw_real X, fftw_real Y,int Mstage) 
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int: x; y, Wl 2, 112; 
long ik; 
fftw_real alpha, lz. wx2, wy2: 


Iz=Lx*Z: 


m12=M/2; 
nl2Z=N/2; 


jik=U: 


for (y=0;y<M;v++)//M 
{ 
Wa (yd) yy 
WV eo o- 
for(x=0;x<N;x++)//N 
{ 
WK = (eh sf ioe) NE XEN 


alpha=PI«lz*(wy2+wx2); 


hglik].re* =aprinv(alpha, Mstage); 
1K oe, 


if( progress update(3)){ 
progress hide(); 
return; 


fftw_real OnuralDecoder:: Magnitude(fftw_complex »*field) 


{ 


int 1: 
fftw_real v.isumV, vl, v2: 


Sui Vv =U: 


for(i=0; i < d_CalcNxd_CalcM; i++){ 
v1=field|i].re: 
v2=field[i].im: 
v=sqrt(vl*v1+v2«v2); 
field[i].re=v: 
field[i].im=0; 
sumV+=v: 


C2 
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return sumV/(d_CalcN * d_CalcM); 


} 


void OnuralDecoder::uchar 2field() 


{ 


int arraysize=d_CalcN*d_CalcM*sizeof(fftw_complex ); 
int fieldpos, ucharpos, x, y: 

uchar minval, maxval: 

fitw_real scale; 


d_Field=(fftw-complex *)calloc(arraysize. 1);//calloc clears to zero... 
minval=0; 
maxval=259; 


for(y=0; y< dImageM; y++ ){ 
for(x=0; x<d_ImageN; x++){ 
ucharpos=3*(x + yxd_ImageM); 
minval=MIN(minval, d_data[ucharpos]): 
maxval=MAX(maxval, d_data[ucharpos]): 


j 
} 
scale=1.0/(fftw_real)(maxval-minval),; 


for( y=0; y< d_tImageM; y++ ){ 
for( x=0; x<d_ImageN; x++){ 


hieldpos=x + y*d_CalcN; 
ucharpos=3*(x + y*d_ImageN); 


d_Field|fieldpos].re=((fftw_real)(d_data[ucharpos]-minval) )*scale: 


} 


j 


uchar *OnuralDecoder::field2uchar() 


{ 


int arraysize=3*d_lmageN*d_ImageM«sizeof(uchar): 
int fieldpos, ucharpos, x, y; 

fftw_real ininval, maxval, v, vre, vim: 

fftw_real scale, val: 


uchar *retval: 


retval=(uchar* )calloc(arraysize, 1);//calloc clears to zero... 
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minval= 1e20; 


maxval=0; 


for(v=0; y< dtImageM; y++ ){ 
for(x=0; x<d_ImageN; x++){ 


} 


fieldpos=x + y*d_CalcM: 
vre=d_Field|fieldpos].re: 
vim=d_Field[fieldpos].im; 


v=sqrt(vrexvre+vim*viln ); 


minval=MIN(minval, v); 
maxval=MAX(maxval, v); 
d_Field[fieldpos].re=v; 


scale=255.0/(maxval-minval); 


for(y=0; y< d_tImageM; y++ ){ 
for(x=0; x<d_ImageN; x++){ 


} 
j 


heldpos=x + y*d_CalcN; 
ucharpos=3*(x + y*d_ImageN); 


val=(uchar) ((d_Field|fieldpos].re-minval) scale); 
retval{ucharpos]=(uchar)val; 
retval{ucharpos+1]=(uchar )val; 

ret val{ucharpos+2]=(uchar)val; 


return retval: 


uchar* OnuralDecoder::mainLoop() 


{ 


fitw_real DC: 


int progresscount=0: 


progress—set Max(5); 
progress—update(0); 
progress reset (); 
progress—show(); 


if( progress—setStatus("Transforming Input...")){ 
progress hide(): 


74 





A.2. DECODER: ONURAL RECONSTRUCTION ALGORITHM. 75 


return 0: 


fttwnd_plan ForwardFFTPlan, BackwardFFTPlan; 
uchar2field(): 


progress update(+-+progresscount ); 
if(progress—setStatus("Making FFTW Plans...")){ 

progress hide(): 

return (): 


ForwardFFTPlan=fftw2d_create_plan(d_CalcN, d_CalcM, 
FFTW_FORWARD, 
FFTW_ESTIMATE|FFTW_IN_PLACE): 


BackwardFFTPlan=fftw2d_create_plan(d_CalcN, d-CalcM. 
FFTW_BACKWARD, 
FFTW_ESTIMATE|FFTW_IN_PLACE), 


DC=4«Magnitude(d_Field); //The 4 handles the doubling of arraysize 


noDClquad(d_Field, DC); 


progress— update(+-+ progresscount ); 
if(progress—setStatus("Calculating Forward transform...")){ 
progress— hide(); 

return 0: 


fftwnd_one(ForwardFFTPlan, d_Field, 0); 
d_Field{0].re=0; 
d_Field{0].im=0; 


progress— update(+-+progresscount ); 
if(progress—setStatus("Starting Onural reconstruction...")){ 

progress— hide(); 

return 0: 


ie 


onural(d_Field, d-CaleN, d-CalcM, d_ImageLambda, d_decodeDist, 
dtImageW, d_ImageH, 5); 
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progress update(+-+ progresscount ); 
if(progress—setStatus("Calculating the Reverse Transform...")){ 

progress hide(): 

return 0): 


fftwnd_one(BackwardFFTPlan, d_Field, 0); 


progress update(++progresscount ); 
if( progress—setStatus("Complete")){ 

progress hide(); 

return (0); 


\. 
uchar+ decoded=feld2uchar(); 
free(d_Field); 


progress hide(); 
return decoded; 
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